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Abstract 

Using the density-matrix renormalization-group method we study 
the surface critical behaviour of the magnetization in Ising strips in 
the subcritical region. Our results support the prediction that the 
surface magnetization in the two phases along the pseudo-coexistence 
curve also behaves as for the ordinary transition below the wetting 
temperature for the finite value of the surface field. 
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1 Introduction 



More than two decades of recent studies have yielded a fairly detailed un- 
derstanding of the critical behavior at surfaces ^1 12]. However, attempts 
to verify theoretical predictions, both in experiments and in model systems, 
often point out issues which need further clarifications. We came across such 
an issue in the recent work by Brovchenko et al., where the surface critical 
behavior of a model water in the slitlike and cylindrical pores |3,, and of a 
Lennard- Jones fluid in the slitlike pores j3] was studied by means of Monte 
Carlo simulations. In both cases fluid particles were assumed to interact with 
a wall via a (10-4) long-range potential and a parameter that measures the 
well depth of the wall-fluid potential was chosen to correspond to a weakly 
attractive surface. A one-component fluid like water or the Lennard- Jones 
fluid is expected to lie in the universality class of the normal surface transition 
of semi-infinite Ising system. In a magnetic language the normal transition 
is characterized by two relevant scaling fields, a surface scaling field c > 
and a non-zero external surface field \hi\ > 0. c describes the enhancement 
of interactions in the surface layer, c > corresponds to a reduced tendency 
to order in the surface, which is the case generic for fluid systems because 
the presence of a wall should decrease the net fluid-fluid attraction between 
a molecule and its nearest neighbors below the corresponding bulk value. On 
the other hand, the containing walls exert an effective potential on a fluid 
and in magnetic language this coresponds to some generally nonzero surface 
field hi. There is a possibility to mimic the situation of vanishing surface field 
hi = 0, i.e., the so called ordinary surface transition behavior, by suitable 
tunning wall-fluid interactions relative to fluid-fluid interactions. Due to the 
lack of Ising symmetry in a "real" fluid, it is very unlikely to find a wall-fluid 
potential that corresponds exactly to /ii = 0, however one does find a wall- 
fluid potential which is "neutral". As was shown in Ref. JE\ for T > Tc, this 
"neutral" wall gives rise to the Gibbs adsorption F ~ that is constant along 
the critical isochore and is characterized by a fluid density profile which, away 
from the walls where oscillations arise, is almost flat throughout the slit [S]. 
For the "neutral" wall a parameter that measures the well depth of the wall- 
fluid potential corresponds to a weakly attractive surface. In Refs. [21111 even 
more weakly attractive substrates were considered. The authors focused on 
the subcritical regime and studied the temperature dependence of density 
profiles along the pore liquid-vapour coexistence curve. Recall, that the nor- 
mal transition is governed by the fixed point of the renormalization-group 
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transformation hi = oo,c = oo and should be equivalent to the extraordinary 
transition given by fixed point hi = 0,c = —oo jEllZllHl- At these (equivalent) 
surface transitions, the order parameter (OP) at the surface layer mi should 
have a leading thermal singularity of the same form as the bulk free energy, 
i.e., \Tc — Tp~°, where Tc is bulk critical temperature. More precisely one 
expects for r ^ [7| |H] the limiting behavior 

mi - {mic + AiT + A2T^ + . . .) ^ (1) 

where r = (Tc — T)/Tc, and the contribution in parentheses is a regular 
background. For both model fiuids Brovchenko et al. [SI 1^ defined the local 
OP as Ap{z) = {pi{z) — pv{z))/2, where pi{z) and Pv{z) are the density 
profiles of the coexisting liquid and vapour phases, respectively, and found 
that below the bulk critical temperature Tc this OP shows the behavior which 
is in accordance with the ordinary transition. In particular, near the surface 
a variation of Ap with reduced temperature r follows the scaling law with a 
value of the exponent close to the /3i ~ 0.82 of the ordinary transition in the 
Ising system in d=3, i.e. 

Api(r)^r^i. (2) 

On the basis of these observations, made for the confined fluids, the authors 
put forward the hypothesis that the difference Ap between the densities of 
coexisting phases near the surface should follow the behavior (j21) also near 
strongly attractive surfaces below the wetting temperature Tyj. This is based 
on the assumption that the term ~ r^^ should always be present in both coex- 
isting phases below T^. The authors reconsider the surface critical behavior 
of the semi-infinite Ising model claiming that below the wetting tempera- 
ture the surface magnetizations m{ and in the two phases along the 
coexistence curve should have the following limiting behavior for r — >^ 0: 

mi = BiT^' + mic + A'lT + . . . + A^^JtI^-'' (3) 
mf = -Birf^' + mic + A[t + . . . + A^_Jt\^-'' (4) 

The symmetric term ~ which describes the temperature dependence 
of the magnetization at /ii = 0, accounts for the missing-neighbor effect 
and, as the authors claim, was overlooked in Ref. [71|H]. Above the wetting 
temperature there exist a single phase which is expected to have the surface 
magnetization of the form given by the above equation. 

For d = 2 semi-infinite Ising model there exist exact results for mi in 
the presence of the surface field. They were derived by McCoy and Wu and 
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also by Au-Yang and Fisher using a Pffafian method GDI- Specifically, 
let us consider a planar rectangular lattice with coordinates i (horizontal) 
and j (vertical) with spins aij = ±1 located at the sites of the lattice and 
interacting with nearest-neighbors via the coupling K = f3J > 0, = l/fc^T. 
Assume vanishing bulk magnetic field /i = 0, a cycling boundary condition 
in the horizontal direction and a surface magnetic field hi, measured in the 
units of the coupling constant J, interacting with one of the two horizontal 
boundary rows of spins. On the second boundary the spins are free. The 
configurational Hamiltonian is defined as 



(5) 



where the sums run over 1 < i < M and 1 < J < A^. The analytic expression 
for the free energy of this system was obtain as a sum |5] 



MNF + 2NTo + NJ^. 



(6) 



The first term is the bulk free energy, and the terms 2NJ-'o and NJ-' are ad- 
ditional contributions coming from the existence of the free boundary which 
interacts with the surface magnetic field. The magnetization mi of the first 
row was calculated from 

= -dhi 

and various limiting cases were discussed. In the case relevant for the ordi- 
nary transition, i.e., for hi —>■ 0, the boundary spontaneous magnetization 
exists only in the thermodynamic limit M, N oo: 



mi(0^ 



lim mi{hi) 



cosh 2/5 J-coth 2(3 J 



cosh 2/? J 

This vanishes at the critical temperature as 

"21n(l + v^)^'^' 



1/2 



(8) 



mi(0^ 



72-1 



1/2 



(9) 



from which one can read off the value of the surface critical exponent P"'^'^ 
for d = 2 Ising model, i.e., Z?"^"^ = 1/2. The case when T is near Tc and hi 
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is positive and away from zero is relevant for the normal transition. For this 
case the result is: 



mi{h\ 



Taylor ser. in r 




In Irl 



(10) 



where zi = tanh/5/ii. Thus the leading singularity of the boundary mag- 
netization agrees with the prediction by Diehl for the normal surface 
transition, i.e., mi has a leading thermal singularity of the same form r^"" 
as the bulk free energy, a = in ci = 2 Ising model which corresponds to the 
logarithmic behavior. In Ref. the first two terms of the Taylor series in 
r were given explicitly 



where 



2K r 



(11) 



:TxDAzA 



+ (1 + V^)' 



l + {l + V2fzl\nzl 
1 



-77 



+ y2 + ln(l + l/V2 



(12) 



On the other hand, if hi is nonzero but sma// then, for T near Tc, the limiting 
behavior of mi is different: 



mi 



^1/2 



(13) 



where a = (1 — z)/z{l + z), z = tanh/T and a = 1 —AK^t + 0{t'^) for r ^ 0. 
Thus, as Zl Eq. (fT^ agrees with Eq. (jH)) and exhibits the square-root 
behavior of the ordinary transition. These exact results show explicitely 
that for strong surface field the prediction Q is not true sufficiently close 
to the critical temperature. On the other hand in view of Eq. (fTSj) it is 
understandable why simulation results for weakly attracting substrates [SI Ej 
may show the ordinary transition behavior. However, results described above 
concern the behavior of the boundary magnetization only for one of the 
two possible bulk phases, and the temperature dependence of the difference 
between the magnetizations of both coexisting phases near the surface has not 
been studied. This is due to the fact that in the absence of the bulk magnetic 
field the choice of the sign of hi breaks the symmetry in the finite system. 
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and, for example, the positive surface field yields (+) phase in the bulk in 
the thermodynamic limit. In order to calculate the boundary magnetization 
for the case of the (— ) bulk phase in the presence of the positive boundary 
field, one would have to perform calculations in the presence of infinitesimaly 
small negative bulk field and put h —>■ 0~ after taking the thermodynamic 
limit or to solve the model with very sophisticated boundary conditions. 



A 




Figure 1: Schematic phase diagram for an Ising strip with positive surface 
fields. There are two thermodynamic paths presented: a) along the bulk 
coexistence and b) along the (pseudo) coexistence of the confined system. 

So far exact solution of the d = 2 Ising model at nonvanishing bulk 
magnetic field is not available, however, recently developed density-matrix 
renormalizat ion-group (DMRG) method ^1] allows for very accurate numer- 
ical calculations in the presence of the arbitrary surface and bulk fields. The 
DMRG method, based on the transfer matrix approach for calculating the 
partition functions, includes the critical fluctuations and therefore is very 
suitable for studying the critical behavior. In the following we will use this 
approach to calculate the magnetization profiles in the strip geometry, i.e., 
the geometry analogous to the one studied in Refs. O Ej . We assume iden- 
tical surface fields /ii = /i2 > acting on the two boundaries separated by a 
finite distance L and consider two different thermodynamic paths: (i) along 
the bulk coexistence h = 0, and (ii) along the pseudo-coexistence of the 
confined system h = hco{T) (see Fig^. The first path is the one for which 
exact results summarized above have been obtained. In this case we want 
to explore how the finite-size effects in the confined geometry influence the 
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crossover from the one type of the surface critical behavior to another. The 
second path allows to study the temperature dependence of the difference 
between the surface magnetizations of both pseudo-coexisting phases. 

The DMRG provides an efficient algorithm to construct the effective 
transfer matrix Tc for large two-dimensional classical systems at finite tem- 
peratures ^2]- Starting with a small system (e.g. L = 4 in our case), for 
which Tc can be diagonalized exactly, one adds iteratively couples of spin 
columns until the allowed (in the computational sense) size of the effective 
matrices is reached. Then further addition of new spins forces one to discard 
simultaneously the least important states to keep the size of the effective 
transfer matrices fixed. This truncation is done through the construction of 
a reduced density matrix whose eigenstates provide the optimal basis set m^. 
The size of the effective transfer matrix is then substantially smaller than the 
original dimensionality of the configurational space (2mA) ^ 2^. Generally, 
the larger is mx, the better accuracy is guaranteed. In the present case, we 
keep this parameter up to mx = 40. Typically a truncation error was not 
larger than 10~^^. We estimate that the errors in the plots are smaller than 
the symbol size. The DMRG method allows to study much larger systems 
(up to L = 340 in this paper) than it is possible with standard exact di- 
agonalization method which can handle with systems up to several dozens 
columns for Ising model. Comparisons with exact results for the case of 
vanishing bulk magnetic field show that this technique gives very accurate 
results in a wide range of temperatures [13]. 

2 Along the bulk coexistence h = 

First we discuss results for h = 0. Recall, that in a finite system with the 
positive hi and h = 0, below Tc there exist only a (+) phase characterized 
by magnetization profiles m\z) which are positive across the strip [14 . In 
the Fig. we show the log-log plot of the difference m{ — as a function 
of r calculated for the strip of the width L = 340 and for the selection of 
the surface fields. For the weakest considered surface fields, i.e., for hi = 
0.001, 0.005, 0.01, and 0.03, we find the square-root behavior of the ordinary 
transition but only for temperatures r > 0.01. In an agreement with the 
exact result (fTSj) the amplitude of this leading decay does not depend on 
hi. For smaller r there is a crossover to the linear behavior with the hi- 
dependent amplitude. The range of temperatures in which the crossover 
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Figure 2: The difference (m( - m{^) as a function of r at L = 340 for 
various surface fields (circles /;,i=0.001, squares /ii=0.005, diamonds /ii=0.01, 
triangles /ii=0.03, crosses hi=0.1, stars /ii=0.5): a) the log-log plot b) the 
effective exponent. The dashed line denotes the slope 1/2 and the dotted- 
dashed line describes the slope 1. 

takes place depends sensitively on the value of the surface field, the weaker 
hi the further away from the crossover starts, but in any case the linear 
behavior is observed for r < 0.001. This is very well illustrated in the plot of 
the effective exponent of m{ — as a function of r (Fig. |2l3)- The effective 
exponent, i.e. the following quantity 



is the discrete derivative of the data in the log- log scale plot. Such quantity 
probes the local slope (at a given reduced temperature T(i)) providing a 
better estimate of the leading exponent than a log- log plot. The calculation 
of Zi requires very accurate data that can be quaranteed by DMRG data. 

The crossover region is associated with the formation of the maximum 
of the local exponent; it shrinks as the surface field becomes stronger and 
disappears altogether for hi = 0.1. For the strongest considered hi, i.e., 
for hi = 0.5 we find a linear behavior for r < 0.005; again the amplitude 
of this leading decay depends on the surface field. Our findings are consis- 
tent with the exact results (fTUI) (fT^ . hi ^ 0.1 being the approximate value 



lnmi(i -|- 1) — Inmi(z) 



(14) 



lnr(i -|- 1) — Inr(z) 
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Figure 3: The difference (m( - m{^) as a function of r at /ii= 0.001 for various 
strip widths L (circles L =50, squares L =100, diamonds L = 200, triangles 
L = 340): a) the log-log plot b) the effective exponent. The dotted-dashed 
line denotes the slope 1, whereas the dashed describes the slope 1/2. 



of the surface field for which one type of the limiting behavior crosses over 
to another. We are not able to decide whether the crossover to the linear 
behavior observed for weak surface field is connected with entering the su- 
percritical region above the pseudocritical temperature Tc^l- Recall, that the 
shift of the critical point due to the finite wall separation and the symmetry- 
breaking boundary conditions is given by (ignoring nonuniversal metric fac- 
tors) ATc = Tc,L-Tc L-^/''XcihiL^'/''),Ahc L-^/^YdhL^^/") with 

z/ = 1, A = 1/15 and Ai = 1/2 for d = 2 Ising model [Hj. Xc and Yc are 
universal scaling functions. Mean-field analysis show that the scaling func- 
tion Xc is of the order 0(1) for small arguments and very weakly depends 
on the value of hi Thus in mean-field ATc ~ —0.001, but the scaling 
function Xc(C), where C = hiL^^^'^, is not known for d = 2 Ising magnet. 
In principle it can vary strongly with the argument. Although the crossover 
depends sensitively on the value of the surface field, it is not connected with 
the wetting because = (T^ -T^(/ii))/T, ~ 8* 10"^ for = 0.001, ~ 10"^ 
for hi = 0.005, ~ 5 * 10"^ for hi = 0.01, ~ 5 * 10"^ for hi = 0.03. 

Figs. and b show the influence of the finite width of the strip on the 
critical behavior of m( — for the weak surface field, i.e., for hi = 0.001. 
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Figure 4: The magnetization as a function of r at L = 340 for various strip 
layers. The inset presents the enlargement of the large-tau part. 

For larger L the crossover region from the square-root to the linear behavior 
is narrow and located closer to Tc. Notice, that for small systems (see the 
curve for L = 50 in Fig. Eb) the finite-size effects are so strong that the 
ordinary transition behavior is attained only very far from T^. 

In the Fig. 0] we show that the magnetization in the inside layers of the 
wide strip (L = 340) subject to the weak surface fields hi = 0.001 decays as 
~ t'^ with /3 = 1/8, i.e., as the spontaneous magnetization, and then crosses 
over to the linear dependence. Similar behavior for the variation of the local 
order parameter was presented by Brovchenko et al '4^. It is striking that 
the crossover to the linear behavior takes place in approximately the same 
temperature range around ~ 0.01 for all layers. 

3 Along the pseudo-coexistence 

Now we consider the path along the pseudo-coexistence of the confined sys- 
tem h = hco{T] L, hi). When the Ising system is confined between parallel 
walls subject to identical surface fields hi, there is a shift of the bulk first- 
order transition to a finite value of the bulk magnetic field. In order to restore 
the coexistence the sign of the bulk magnetic field h has to be opposite to 
the sign of /ii. In = 2 for the system with finite L there is no unam- 
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biguous way to determine the pseudocoexistence line. One can use several 
criteria, for example, maxima of the specific heat, minima of the inverse cor- 
relation length or inflection points of the solvation force [TS' . However, above 
some characteristic temperature, which corresponds to the (pseudo) capillary 
critical point, the curves based on different criteria separate because they 
are governed by different critical exponent. Here we adapt a very natural 
criterion of the zero total magnetization, i.e., for the fixed value of h one cal- 
culates the total magnetization T = where m/ is the magnetization 
in the /-layer corresponding to a perpendicular distance from the first wall, 
for different temperatures and identifies h = hco{T, hi) at the temperature 
at which F = 0. This method works very good away from the immediate 
neighborhood of the critical point, where the difference between two phases 
vanishes and it is difficult to locate the point corresponding to F = 0. In the 
end we are not able to determine the difference Ami = — too close 
to the bulk critical temperature (in the limit of r ^ 0). 
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E 0.0 
^ -0.5 
-1.0 

1.75 2.00 2.25 
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1.0 
6 0.5 

e" 

0.0 
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T 

Figure 5: The temperature dependence of the surface magnetizations of 
the two coexisting phases calculated along the (pseudo) coexistence lines: a) 
/ii=0.01 b) hi= 0.5. 

We have performed calculations for the strip of the width L = 340 and 
three different surface fields hi = 0.01, 0.5, 0.8. The temperature dependence 
of the surface magnetizations of the two coexisting phases m{ and cal- 
culated along the line h = hco{T, hi) for hi = 0.01 and 0.5 are shown in 
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the Fig. El For the strongest field one can see the asymmetry due to the 
positive surface field. The temperature at which these two curves coincide 
may be identified with Tc^l. Note, that the pseudocoexistence temperature 
is located approximately at the same temperature at which the numerical 
errors become relevant, i.e., for hi = 0.01 at ^ 2.22, for hi = 0.5 at ~ 2.165, 
(Tc ~ 2.26919 for the d = 2 Ising model). The appearence of the numerical 
errors is connected with large fiuctuations near T^^l 




Figure 6: The log-log plots of Ami as a function of r at L=340 for various 
surface fields: a) /ii=0.01 b) hi=0.5 c) hi= 0.8. The vertical dotted lines 
denote wetting temperatures, the dashed line presents the slope 1/2 and the 
dotted-dashed line the slope 1. 

The log-log plots of the difference Ami = m{ — m(^ versus r for hi = 0.01, 
hi = 0.5 and hi = 0.8 are shown in Fig.|nti',b,c, respectively, and the effective 
exponent is presented in Fig. [7| We also mark the wetting temperature 
Tw{hi). The general behavior which can be read off from these plots is that 
below Tyj{hi) the difference Ami = t^i —m{^ decays approximately like t^/^, 
then there is a crossover regime connected with the wetting transition, closer 
to Tc the linear behavior dominates, and finaly there is a rapid decay due 
to the proximity of the pseudo-critical point. For hi = 0.8 the approximate 
square-root behavior takes place in a very narrow range of temperatures, 
because lies quite far away from Tc ((Tc — Tu,)/T^ ^ 0.38). For hi = 0.01 
the linear behavior is not reached since T^(O.Ol) ~ 2.26906. 
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Figure 7: The effective exponent of Ami for various surface fields. 

It is instructive to see the near surface behavior of the magnetization 
profiles in both phases for different temperatures (see Fig. |Ht,b). In the 
phase opposite to the one that is favored by walls, i.e., in the negatively 
magnetized phase m^^, the wetting transition manifests itself by a change in 
dm^^ {l)/dl from the monotonuous function of the distance from the surface 
/ to the one exhibiting a minimum (see Fig. 

In conclusion, our DMRG results for the case of vanishing bulk field 
are in agreement with the exact results and show two different type of the 
asymptotic behavior of the surface magnetization. For weak surface fields 
we see the square-root r— dependence characteristic of the ordinary surface 
universality class which crosses over to the linear behavior sufficiently close 
to Tc. This crossover is not connected with the wetting. For hi > 0.1 we find 
the linear behavior which dominates over the singular r^~", characteristic of 
the normal universality class. Results along hco{T) suggest that Ami ~ t^^^ 
below wetting transition but Ami ~ t above it. However, in order to clear- 
out this issue the exact calculations of Ami in the semiinfinite system is 
needed. Because for very weak surface field the wetting temperature lies 
very close to Tc one may in such a case observe only the square-root behavior 
as in the simulation of the Ref. [30]. 
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Figure 8: The magnetization profiles near the (pseudo) critical point calcu- 
lated along the zero- magnetization line: a) b) m^^ . 

References 

[1] For a general review of critical behavior at surfaces, see K. Binder, 
Phase Transitions and Critical Phenomena edited by C. Domb and 
J. L. Lebowitz (Academic Press, London, 1983), Vol. 8, p. 1. 

[2] H. W. Diehl, Phase Transitions and Critical Phenomena edited by 
C. Domb and J. L. Lebowitz (Academic Press, London, 1986), Vol. 10 
p. 75 and references therein. 

[3] I. Brovchenko, A. Geiger and A. Oleinikova, J. Phys.: Condens. Matter 
16, S5345 (2004). 

[4] I. Brovchenko, A. Geiger and A. Oleinikova, Eur. Phys. J. B 44, 345 
(2005). 

[5] A. Maciolek, R. Evans, and N. B. Wilding, J. Chem. Phys. 119, 8663 
(2003). 

[6] A. J. Bray and M. A. Moore, J. Phys. A 10, 1927 (1977). 

[7] T. W. Burkhard and H. W. Diehl, Phys. Rev. B 90, 3894 (1994). 




14 




Figure 9: The first derivative of the m^^ profiles with respect to the dis- 
tance from the surface at hi=0.5. The strip width is L—3A0. The wetting 
temperature for this surface field is Ty^ ~ 1.96. 

[8] H. W. Diehl, Ber. Bunsenges. Phys. Chem. 98, 466 (1994). 

[9] B. M. McCoy and T. T. Wu, Phys. Rev. 162, 436 (1967). 
[10] H. Au-Yang and M. E. Fisher, Phys. Rev. B 21, 3956 (1980). 
[11] T. Nishino, J. Phys. Soc. Jpn. 64, 3598 (1995). 

[12] Lectures Notes in Physics., edited by 1. Pcschcl, X. Wang, M. Kaulke, 
and K. Hallbcrg (Springer, Berlin, 1999), Vol. 528; U. SchoUwoeck, 
Rev. Mod. Phys. 77, 259 (2005). 

[13] A. Maciolek, A. Ciach, and A. Drzewihski, Phys. Rev. E 60, 2887 (1999). 

[14] M. E. Fisher and H. Nakanishi, J. Chem. Phys. 75, 5857 (1981); ibid. 
78, 3279 (1982). 

[15] A. Maciolek, A. Drzewihski, and R. Evans, Phys. Rev. E 64, 056137 
(2001). 



15 



